package cmp.compile;

import java.sql.ResultSet;
import java.util.HashMap;
import java.util.HashSet;

import cmp.compile.CompileMain;
import cmp.database.*;
import assem.DBConn;
import cmp.main.Globals;
import util.Debug.ErrorReport;

public class Pairwise {
	boolean dolog2 = false; 
	
	public Pairwise(DBConn dbc, String methodPrefix) {
		try {
			long starttime = CompileMain.getTime();
			CompileMain.PrtSpMsg(1, "Update clusters " + methodPrefix + " with PCC results");
			mDB = dbc;
			addPCCtoMethod(methodPrefix);
			CompileMain.PrtSpMsgTime(1, "Complete PCC ", starttime);
		}
		catch (Exception e) {
			ErrorReport.die(e, "Error updating pairwise");
		}
	}
	
	private void addPCCtoMethod(String mPre) {
		try {
			ResultSet rs;
			int [] PGid;		// list of group ids for method
			
	// get library column headers and make SQL and list
			DBinfo info = new DBinfo(mDB, true);
			String [] libList = info.getSeqLib();
			String libSQL = info.getSeqLibSQL();
			
	// get list of POGs from this method
			int PMid = mDB.executeCount("select PMId from pog_method where prefix='" + mPre + "'");	
			int nPOG = mDB.executeCount("select count(*) from pog_groups where PMid=" + PMid);			
			CompileMain.PrtSpMsg(2, "Number of Clusters to process " + nPOG);
			
			PGid = new int [nPOG+1]; // PGid counts 1 to N (not 1 to N-1)
			rs = mDB.executeQuery("select PGid from pog_groups where PMid=" + PMid);
			int x=0;
			while (rs.next()) PGid[x++] = rs.getInt(1);
	
			int goodPCC=0, goodclusters=0;
			long count=0;
	// for each cluster, get its members and libs, compute PCC and save
			for (int k=0; k<nPOG; k++) { // for cluster k
			// get all members		
				HashMap <Integer, String> seqName = new HashMap <Integer, String> ();
				HashMap <Integer, Double []> seqLib = new HashMap <Integer, Double []> ();
				int [] ids;
				int cnt=0;
				
				rs= mDB.executeQuery("select UTid,UTstr from pog_members where PGid=" + PGid[k]);
				while (rs.next()) seqName.put(rs.getInt(1), rs.getString(2));
				
			// get their RPKMs
				ids = new int [seqName.size()]; // for double loop
				double log2 = Math.log(2);
				for (Integer id: seqName.keySet()) {
					rs = mDB.executeQuery("select " + libSQL + " from unitrans "+
							"where UTid=" + id);
					if (!rs.next()) CompileMain.die("Error getting RPKM for " + id);
					Double [] rpkm = new Double [libList.length];
					for (int j=0; j< libList.length; j++)  {
						double r = rs.getDouble(Globals.PRE_LIB + libList[j]);
						if (dolog2 && r != -1) {
							double l = Math.log(r)/log2;
   		    	   				if (l < 0) l = 0;
   		    	   				rpkm[j] = l;
						}
						else rpkm[j] = r;
					}
					seqLib.put(id, rpkm);
					ids[cnt] = id;
					cnt++;
				}
			// compute PCC 
				double maxPCC=-100.0, pcc;
				int cntPCC=0;
	
				for (int i=0; i<cnt-1; i++) {
					int id1 = ids[i];
					for (int j=i+1; j<cnt; j++) {
						count++;
						if (count % 1000 == 0)
							System.err.print(" Processed " + count + " pairs. Processing cluster #" + k + " with " + 
									seqName.size() + " members\r");
						int id2 = ids[j];
						if (id1 > id2) {int t = id1; id1 = id2; id2 = t;}
						
						rs = mDB.executeQuery("select PCC from pairwise " +
								" where UTid1=" + id1 + " and UTid2=" + id2);
						if (rs.next()) pcc = rs.getDouble(1);
						else {
							pcc = getCorrelation(id1, id2, seqLib.get(id1), seqLib.get(id2));
							if (pcc >= Globals.PCCcutoff) {
								mDB.executeUpdate(
										"insert into pairwise (UTid1, UTid2, UTstr1, UTstr2, pcc)" +
										" values (" + id1 + "," + id2 + 
										",'" + seqName.get(id1) + "','" + seqName.get(id2) + "'," + pcc + ")");
							}
						}
						if (pcc >= Globals.PCCcutoff) cntPCC++;
						if (pcc > maxPCC) maxPCC = pcc;				
					}
				}	
				
				int pairs = (cnt * (cnt-1))/2;
				double p = ((double)cntPCC/(double)pairs) * 100.0; 
				if (p>=100) goodclusters++;
				goodPCC += cntPCC;
				
				mDB.executeUpdate("update pog_groups " +
						" set LperPCC=" + Globals.saveDouble(p) + 
						",    LmaxPCC=" + Globals.saveDouble(maxPCC) + 
						" where PGid=" + PGid[k]);
				seqName.clear(); seqLib.clear(); 
			}
			rs.close();
			CompileMain.PrtSpMsg(2, "PCC>=" + Globals.PCCcutoff + " = " + goodPCC  +
					"  Clusters with 100% good PCC = "+goodclusters + "                          ");
		}
		catch (Exception e) {
			ErrorReport.die(e, "Error add PCC to Method");
		}
	}
	
	/**********************************************************
	 * compute all at once but takes a long time
	 */
	public Pairwise(DBConn dbc) {
		try {
			CompileMain.PrtDateMsg("\nAdding pairwise data ");
			mDB = dbc;
			addBlast();
			computePCC();
		}
		catch (Exception e) {
			ErrorReport.die(e, "Error updating pairwise");
		}
	}

	private void addBlast() {
		// read first hit for each pair into database
	}
	/**********************************************************
	 * Compare all pairs. If hasBlast or PCC > 0.5 put in database
	 */
	private void computePCC() {
		try {
			CompileMain.PrtSpMsg(1, "Compute PCC values (first time only)");
			ResultSet rs;
			int nUT = mDB.executeCount("select count(*) from unitrans");
			
		// What pairs have blast results? Key = "UTid1:UTid2" with UTid1 < UTid2
			HashSet <String> hasHit = new HashSet <String> ();
			// TODO - read database once code for reading combined.fasta is in
			
		// get library column headers
			rs = mDB.executeQuery("select allSeqLib from info");
			rs.next();
			String libs = rs.getString(1);
			if (libs==null || libs.equals("")) {
				CompileMain.PrtSpMsg(1, "No libraries in database; skip PCC step");
				return;
			}
			CompileMain.PrtSpMsg(2, "Libraries: " + libs);
			String [] libList = libs.split(" ");
			String libSQL = "";
			for (int i=0; i<libList.length; i++) {
				libList[i] = Globals.PRE_LIB + libList[i];
				libSQL += "," + libList[i];
			}
			
		// Read all UTid, UTstr and L_; 
			HashMap<Integer,Double []> seqLib = new HashMap <Integer,Double []>();
			HashMap<Integer, String> seqName = new HashMap <Integer, String> ();					
			long count=0;
			rs = mDB.executeQuery("Select UTid, UTstr " 
					+ libSQL + " from unitrans where aaLen!=0");
			while (rs.next()) {
				int UTid = rs.getInt(1);
				String UTstr = rs.getString(2);
				Double [] rpkm = new Double [libList.length];
				for (int i=0; i< libList.length; i++)  rpkm[i] = rs.getDouble(libList[i]);			
				seqLib.put(UTid, rpkm);
				seqName.put(UTid, UTstr);
				count++;
			}
			long max = (count/count-1)/2;
			CompileMain.PrtSpMsg(2, "Computing PCC for " + count + " sequences and "  + libList.length + " libraries");
			count=0;
			for (int id1=0; id1<nUT-1; id1++) {
				if (!seqLib.containsKey(id1)) continue;
				for (int id2=id1+1; id2<nUT; id2++) {
					if (!seqLib.containsKey(id2)) continue;
					
					Double [] l1 = seqLib.get(id1);
					Double [] l2 = seqLib.get(id2);
					if (l1 == null || l2 == null) continue;
					count++;
					if (count % 1000 == 0) {
						long xx = (count/max) * 100;
						System.err.print(xx + "%  PCC values done.....\r");
					}
					
					double pcc = getCorrelation(id1, id2, seqLib.get(id1), seqLib.get(id2));
					String key = id1 + ":" + id2;
					if (hasHit.contains(key)) {
						// update
					}
					else if (pcc > 0.5) {
						mDB.executeUpdate("insert into pairwise (UTid1, UTid2, UTstr1, UTstr2, pcc)" +
						" values (" + id1 + "," + id2 + 
							",'" + seqName.get(id1) + "','" + seqName.get(id2) + "'," + pcc + ")");
					}
				}
			}
		}
		catch (Exception e) {
			ErrorReport.reportError(e, "Error computing PCC");
		}
	}
    // Pearson Correlation Coefficient: verified with on-line PCC calculators
	// also in jpave.util.ui.sortable
    private double getCorrelation(int ix, int iy, Double [] xVect, Double [] yVect) {
    		if (xVect==null || yVect==null) {
    			System.err.println(ix + "," + iy + " null xVect or yVect ");
    			return 0;
    		}
    
        double sumX=0.0, sumY=0.0, sumXY = 0.0, powX=0.0, powY=0.0;
        int sz = 0;
        try {
	        for(int i = 0; i < xVect.length; i++)
	        {
	            double x = xVect[i];
	            double y = yVect[i];
	            if (x == -1.0 || y == -1.0) continue;
        			sz++;
	            sumX  += x;
	            sumY  += y;
	            powX  += (x*x);
	            powY  += (y*y);
	            sumXY += (x*y);
	        }	
	        if (sz==0) return 0;
	        double dsz = (double) sz;
	        double num = sumXY - ((sumX*sumY)/dsz);	        
	        double demX = powX - ((sumX * sumX)/dsz);
	        double demY = powY - ((sumY * sumY)/dsz);
	        double dem = Math.sqrt(demX * demY);
	        double r = (num!=0 && dem!=0) ? (num/dem) : 0;
	        return r;
        }
        catch (Exception e) {
        		System.out.println(xVect.length + " :" +  join( xVect) );
        		System.out.println(yVect.length + " :" +  join( yVect) );
        		ErrorReport.reportError(e, "Internal error: computing correlation");
        		return 0.0;
        }
      }//end: GetCorrelation(X,Y)
	
	private String join(Double [] s) {
		
		StringBuffer buffer = new StringBuffer();
		
		for (int i=0; i< s.length; i++) buffer.append(s[i] + " ");
		
		return "'" + buffer.toString() + "'";
	}
	DBConn mDB;
}
